home *** CD-ROM | disk | FTP | other *** search
/ BBS in a Box 3 / BBS in a box - Trilogy III.iso / Files / Prog / U-Z / VideoToolBox Folder / Utilities / CalibrateLuminance / CalibrateLuminance.c next >
Encoding:
Text File  |  1993-03-04  |  37.3 KB  |  972 lines  |  [TEXT/KAHL]
  1. /*
  2. CalibrateLuminance.c
  3. by Denis Pelli, Lan Zhang, and Preeti Verghese
  4. You will need "Numerical Recipes in C" in order to compile this file; see
  5. limitations below.
  6.  
  7. USE
  8.  
  9. You must run this program (or your own equivalent) to calibrate your video card,
  10. ISR Video Attenuator, and monochrome monitor, as described by Pelli & Zhang
  11. (1991).
  12.  
  13. D.G. Pelli and L. Zhang (1991) Accurate control of contrast on microcomputer
  14. displays. Vision Research, 31:1337-1360.
  15.  
  16. The program's results are stored in a file called LuminanceRecord?.h that you
  17. can then use in the rest of your C programs. "?" will be a digit indicating
  18. which monitor. LuminanceRecord?.h describes the gains of the three video
  19. pathways of your video card and ISR Video Attenuator, and describes the gamma
  20. function of your monitor, to allow automatic gamma correction later, using
  21. SetLuminance(), etc.
  22.  
  23. There are two ways of using the LuminanceRecord file. You can use the
  24. preprocessor #include statement, in the midst of your C code, provided you've
  25. got a LuminanceRecord structure called LR. Or you can use the VideoToolbox
  26. routine called ReadLuminanceRecord() to read the data (at runtime) into a
  27. LuminanceRecord of any name.
  28.  
  29. For historical reasons this program not only measures the gamma function, but
  30. also fits a variety of functions to it. In fact the latest implementation of
  31. Luminance.c, which is the only routine that makes direct use of the gamma
  32. description, would be quite content with just the tabulated gamma function. At
  33. present the gamma function is tabulated very coarsely by CalibrateLuminance, so
  34. it doesn't bother to put this table into the include file LuminanceRecord.h.
  35. However, if you'd rather write your own substitute for this program, then I
  36. would suggest just measuring the luminance at 256 different levels (from 0 to
  37. VMax) and saving the table in the include file. You will also need to measure
  38. the RGB gains, but that can easily be done directly, using an oscilloscope,
  39. rather than trying to be clever, as in here, where we invert the nonlinearity to
  40. infer voltages from the nonlinear luminances that we measure.
  41.  
  42. CalibrateLuminance also produces a CricketGraph file, CalibrateLuminance%d.data
  43. (where %d is replaced by the monitor number), which is suitable for graphing the
  44. monitor's gamma function. Use the CricketGraph format file
  45. CalibrateLuminance.format.
  46.  
  47. This program measures the luminance of a patch on the screen, using a photometer
  48. and, optionally, a 12-bit Analog to Digital Converter (ADC). If you have a Data
  49. Translation Forerunner Analog-to-Digital card, then GetVoltage.c will
  50. automatically find and use it. You can set the MANUAL flag to force manual
  51. operation. If your photometer has an analog output, it's convenient and more
  52. accurate to have all the voltages read in automatically via an ADC, but you
  53. won't need to run CalibrateLuminance very often, and you can get by with manual
  54. calibration.
  55.  
  56. You should set the background luminance to approximately the same value as you
  57. will use in your experiment. Note that the channels gains depend on the DACs as
  58. well as the ISR Video Attenuator, and will vary from DAC to DAC by ±5%. That's
  59. why you must do this calibration yourself using your own video card and monitor.
  60.  
  61. This program can calibrate any of your screens, including the main screen. It
  62. now works fine with a single screen, alternating the user dialog with the
  63. measurements.
  64.  
  65. LIMITATIONS
  66.  
  67. CalibrateLuminance uses routines from Numerical Recipes in C to do the
  68. polynomial and power law fits to the gamma function. They're copyrighted, so I
  69. can't distribute them. Note: I HAVE included a compiled application
  70. CalibrateLuminance that you can use NOW. You only need to buy the Numerical
  71. Recipes if you want to MODIFY CalibrateLuminance.c. See "Improve Numerical Recipes" 
  72. in the notes folder.
  73.  
  74. HISTORY
  75.  
  76. 4/25/89    Preeti and Denis
  77. 6/18/89    Denis added numerical display of clut index and triple grating for first screen.
  78. 8/4/89    Denis replaced SetEntries call by GDSetEntries, and generally updated everything.
  79. 8/18/89    Lan Measure the whole routine(R+G+B and R G B gain) 40 times to eliminate the
  80.         effect of screen luminance drift.
  81. 9/8/89    Lan & Denis generally tidied it up.
  82. 9/8/89    denis: replaced polynomial fit by powerlaw fit
  83. 9/10/89 denis: got gain measurement to work with sufficient accuracy. 
  84.         This involved many small changes. I now
  85.         wait for a second after any large luminance change to allow the
  86.         photometer to settle. I also measure the channel gains at many different
  87.         settings of the other DACs, in order to average out the effect of DAC
  88.         inaccuracies.
  89. 10/30/89 Lan & Denis: introduced the option of manual readings, and made console smaller.
  90. 11/17/89 Lan & Denis: cleaned up for general release.
  91. 11/30/89 Denis: added comment to LuminanceRecord.h explaining power[] parameters.
  92. 3/29/90    dgp    2.15 Updated to use new GetVoltage that looks for ForeRunner card,
  93.             and new GDOpenWindow() that uses CWindowPtr instead of WindowPtr. Introduced
  94.             conditionals so it compiles without errors under MPW C 3.1. 
  95. 4/21/90    dgp    2.20. Fixed the bug in the fixed-power power law fit. Updated whole file
  96.             to be compatible with latest versions of all subroutines. Corrected printout
  97.             of nBackground.
  98. 7/28/90    dgp    2.3. Added code (in GetALuminance.c) to find the equivalent number to
  99.             produce any desired background luminance.
  100.             All luminance measurement is now done via the new subroutine GetALuminance. 
  101.             Automatic and manual measurements now use the exactly the same code.
  102.             Calibration of (linear) RGB gains is now optional.
  103. 9/18/90    dgp    2.4. Changed all instances of "v" to "V". Pelli
  104.             & Zhang (1991) refer to a nominal voltage v; this file
  105.             now refers to the "equivalent number" V; they are related by V=VMax*v.
  106. 9/22/90    dgp    2.5. Fixed cosmetic errors that prevented use of this program on the main
  107.             screen. The trick is appropriate use of BringToFront() and SendBehind()
  108.             referring only to my own window. Trying to do BringToFront() on the console
  109.             often caused a bus error, for no obvious reason. I'll have to ask Mike
  110.             Kahl.
  111. 9/24/90    dgp    2.6. Added screen to the LuminanceRecord.h file and appended it to
  112.             the LuminanceRecord.h and .data filenames. This makes it easy to calibrate
  113.             all your monitors and keep the records straight.
  114.             Added LR.date string.
  115.             New default is to retain old gainAccuracy when rgb gains are not remeasured.
  116. 10/10/90 dgp Added SetDepth(). Fixed bug that reducing the number of frames sampled
  117.             per a/d measurement by the number of cycles. The number of frames per
  118.             measurement is now fixed.
  119. 10/12/90 dgp Added LR.dpi and LR.Hz to the LuminanceRecord?.h file.
  120. 10/17/90 dgp Removed unused variables. Added reference to paper to LuminanceRecord?.h.
  121. 12/12/90 dgp Measure and subtract off the dark voltage.
  122.             Fixed dumb error of assuming wrong type when printing LP->VMax
  123.             & LP->coefficients, which was resulting in zeroes in the LuminanceRecord.
  124. 4/15/91    dgp Check for NewPaletteManager().
  125. 8/24/91    dgp    made compatible with THINK C 5. This required introducing several
  126.             casts: (unsigned long *) and (unsigned char *). Hopefully, this won't
  127.             compromise compatibility with old THINK C 4.
  128. 3/10/92    dgp    include mc68881.h
  129. 6/23/92    dgp    added quick check of photometer gain.
  130. 8/27/92    dgp    replace SysEnvirons() by Gestalt()
  131. 10/23/92 dgp read latest LuminanceRecord for default values. Mentioned ReadLuminanceRecord()
  132.             in the LuminanceRecord header file.
  133. 10/24/92 dgp use GDFrameRate() to measure the frame rate.
  134. 11/13/92 dgp advertise ReadLuminanceRecord.c
  135. 12/17/92 dgp Enhanced to support arbitrary dacSize. Get dacSize by calling GDGetGamma(). 
  136.             Mostly I just replaced 255 by LP->VMax. Didn't test it after changes.
  137. 12/21/92 dgp No longer load unused dac bits.
  138. 2/23/93    dgp    Use new GDOpenWindow1 and GDDisposeWindow1.
  139. */
  140. #include "VideoToolbox.h"
  141. #include <math.h>
  142. #include "Luminance.h"
  143. #include <Fonts.h>
  144. #include <Packages.h>
  145. #if THINK_C
  146.     #include <console.h>
  147. #endif
  148. #include <nr.h>            /* prototypes of Numerical Recipes and definition of FLOAT */
  149.  
  150. #define AUTOMATIC        1    /*     =0:input reading by hand from the photometer to the terminal
  151.                                    =1: use A/D converter to read from the photometer directly */
  152. #define NLRGB            16    /* number of times to measure each RGB gain, should be >1*/
  153. #define LUMINANCES        32    /* approx. number of different luminances to measure */
  154. #define    NFRAME            60    /* number of frames per ADC luminance measurement */
  155. #define MAX_WINDOWS        6    /* maximum number of screens that you might have */
  156. #define ROUND_LUMINANCES (2+(256+256/LUMINANCES-1)/(256/LUMINANCES))
  157.                             /* round LUMINANCES up to the next divisor of 256, add 2 */
  158.  
  159. void main(void);
  160. void CalibrateLuminance(void);
  161. void SaveData(luminanceRecord *LP,int nL,int V[],double L[]
  162.     ,int nLrgb,int nrgb[3][NLRGB][3],double Lrgb[3][NLRGB],double gain[3]);
  163. void SaveLuminanceRecord(luminanceRecord *LP);
  164. FLOAT PowerRMSError(FLOAT p[]);
  165. double GetALuminance(luminanceRecord *LP,GDHandle device,
  166.     int frames,double LuminancePerVoltage,int entry,int red,int green,int blue);
  167. void LToVHunt(luminanceRecord *LP,GDHandle device,CWindowPtr window,
  168.     double LuminancePerVoltage,int frames,double darkLuminance);
  169.  
  170. /* These variables are out here, as globals, so that PowerRMSError can access them directly */
  171. static int nL=ROUND_LUMINANCES;
  172. static double L[ROUND_LUMINANCES];    /* luminance at V[i] */
  173. static int V[ROUND_LUMINANCES];
  174. static int variables=4;
  175. static FLOAT *p;
  176.  
  177. void CalibrateLuminance(void);
  178.  
  179. void main(void)
  180. {
  181.     Require(gestalt8BitQD);
  182.     CalibrateLuminance();
  183. }
  184.  
  185. void CalibrateLuminance(void)
  186. {
  187.     int error;
  188.     long value;
  189.     int i,j,k,ii;
  190.     static char string[100];
  191.     static Rect myRect,TestRect,labelRect;
  192.     int bottom;
  193.     GDHandle device = NULL,oldDevice = NULL;
  194.     WindowPtr window=NULL,oldWindow=NULL;
  195.     static WindowPtr windows[MAX_WINDOWS];
  196.     short int FontNum;
  197.     int nLrgb=NLRGB;
  198.     int screen;
  199.     static double Lrgb[3][NLRGB];
  200.     static int nrgb[3][NLRGB][3];
  201.     static double gain[3];
  202.     static double LuminancePerVoltage=1000.0;
  203.     static luminanceRecord LR,*LP;
  204.     static double *x,*y,*sig;                        /* for polynomial fit */
  205.     static FLOAT *a,**u,**v,*w,**cvm,chisq;            /* for polynomial fit */
  206.     int ma;                                            /* for polynomial fit */
  207.     static FLOAT **xi,ftol,fret;                    /* for power law fit */
  208.     int iter;                                        /* for power law fit */
  209.     static double e,f,VV;
  210.     int ctSize;
  211.     int readNumber,readTotal;
  212.     Boolean flag;
  213.     short int testSize,cycles,textSize;
  214.     Boolean automatic,skipRGB;
  215.     char myChar;
  216.     int frames;
  217.     double luminance;
  218.     static double darkLuminance;
  219.     GammaTbl *gammaTblPtr=NULL;
  220.     
  221.     automatic=AUTOMATIC && (CardSlot(".ForeRunner")>0);    /* disable if card is missing */
  222.     LP=&LR;
  223.     #include "LuminanceRecord1.h"                         /* get defaults from old calibration */
  224.     #if THINK_C
  225.         console_options.title="\pCalibrateLuminance";    /* change the console window */
  226.         console_options.nrows=12;
  227.         console_options.top=20;
  228.         console_options.left=0;
  229.         console_options.ncols=105;
  230.         console_options.txSize=9;
  231.         printf("\n");                                    /* ask THINK C to initialize QuickDraw */
  232.         CopyQuickDrawGlobals();                            // Make sure qd is valid.
  233.     #else
  234.         InitGraf(&qd.thePort);
  235.         InitWindows();
  236.         InitFonts();
  237.     #endif
  238.  
  239.     oldDevice = GetGDevice();
  240.     GetPort(&oldWindow);
  241.     GetFNum((StringPtr) "\pChicago",&FontNum);
  242.     for(i=0;i<MAX_WINDOWS;i++){
  243.         device = GetScreenDevice(i);
  244.         if(device == NULL)break;
  245.         windows[i] = GDOpenWindow1(device);
  246.         SendBehind(windows[i],NULL);
  247.         SetPort(windows[i]);
  248.         TextFont(FontNum);
  249.         textSize=128;
  250.         TextSize(textSize);
  251.         myRect = windows[i]->portRect;
  252.         bottom=myRect.bottom;
  253.         MoveTo(myRect.right/2,myRect.bottom/2);
  254.         sprintf(string,"%1d",i);
  255.         CtoPstr(string);
  256.         Move(-StringWidth((void *)string)/2,textSize/2);
  257.         DrawString((StringPtr) string);
  258.     }
  259.     SetPort(oldWindow);
  260.     SetGDevice(oldDevice);
  261.  
  262.     printf("Welcome to CalibrateLuminance (version 12/17/92).\n\n"
  263.     "This program helps you to calibrate your ISR Video Attenuator, as well as your video card\n"
  264.     "and monitor. 'Enter' indicates that you should type in a number or letter followed by 'RETURN'.\n"
  265.     "Most questions provide a default answer (in parentheses).\n"
  266.     "Just hit 'RETURN' if you want the default. (Or hit Command-. to quit.)\n"
  267.     "Where appropriate, the default is taken from latest LuminanceRecord.\n");
  268.  
  269.     screen=0;
  270.     if(GetScreenDevice(LP->screen) != NULL)    screen=LP->screen;
  271.     printf("\nEnter the number of the screen you want to calibrate (%d):",screen);
  272.     gets(string);
  273.     sscanf(string,"%d",&screen);
  274.     printf("%d\n",screen);
  275.     for(i=0;i<MAX_WINDOWS;i++){
  276.         device = GetScreenDevice(i);
  277.         if(device == NULL)break;
  278.         GDDisposeWindow1(windows[i]);
  279.     }
  280.     LP->screen=screen;
  281.     sprintf(string,"LuminanceRecord%d.h",screen);
  282.     i=ReadLuminanceRecord(string,LP,0);        // try to read latest luminanceRecord
  283.     if(i<1)printf("Couldn't find an old copy of “%s”.\n",string);
  284.     else printf("Read old “%s” to set default values.\n",string);
  285.         
  286.     GetPort(&oldWindow);
  287.     oldDevice = GetGDevice();
  288.     device = GetScreenDevice(screen);            /* choose which screen to test */
  289.     ctSize=(*(**(**device).gdPMap).pmTable)->ctSize;        
  290.     if(ctSize<3 && NewPaletteManager())
  291.         SetDepth(device,8,1,1);            /* 8-bit pixelSize & color mode */
  292.     ctSize=(*(**(**device).gdPMap).pmTable)->ctSize;        
  293.  
  294.     // Get dacSize
  295.     error=GDGetGamma(device,&gammaTblPtr);        // Takes 200 µs.
  296.     if(error) LP->dacSize=8;
  297.     else LP->dacSize=gammaTblPtr->gDataWidth;
  298.     printf("The video card has %d-bit dacs.\n",LP->dacSize);
  299.  
  300.     window = GDOpenWindow1(device);
  301.     SendBehind(window,NULL);
  302.     printf("Test screen is set to %d colors ",ctSize+1);
  303.     if(ctSize <3)     {
  304.         PrintfExit("\n\007 Sorry. You need to set the test screen to at least 4 colors.\n"
  305.             "Hit 'RETURN' to quit the program then go to 'Monitors' in the "
  306.             "Control Panel and set to 'Colors' and '256'\n");
  307.     }
  308.     GDGetGray(device,&flag);        /* flag=0:color mode;flag=1:gray mode */
  309.     /* request color mode */
  310.     if (flag != 0 && NewPaletteManager())
  311.         SetDepth(device,(**(**device).gdPMap).pixelSize,1,1);
  312.     GDGetGray(device,&flag);        /* flag=0:color mode;flag=1:gray mode */
  313.     if (flag != 0)    {
  314.         PrintfExit("\n\007 Sorry. You need to set the test screen to color mode.\n"
  315.             "Hit 'RETURN' to quit the program. Go to 'Monitors' in the "
  316.             "'Control Panel' and set to 'Colors' and '256'.\n");
  317.     }
  318.     printf("and color mode\n");
  319.     LP->Hz=GDFrameRate(device);
  320.     printf("The test screen's frame rate is %.2f Hz\n",LP->Hz);
  321.     
  322.     myChar='N';
  323.     printf("\nThere are two kinds of measurement:\n"
  324.     "1. Compulsory calibration of the (nonlinear) monitor's gamma function.\n"
  325.     "This will be affected by your choice of background luminance.\n"
  326.     "2. Optional calibration of the (linear) ISR Video Attenuator 3 channel gains.\n"
  327.     "This is independent of the background. It takes about 15 minutes.\n\n"
  328.     "Do you wish to skip 2., the optional three-channel calibration? (%c):",myChar);
  329.     gets(string);
  330.     sscanf((char*)string,"%c",&myChar);
  331.     skipRGB=(toupper(myChar)=='Y');
  332.     if(skipRGB) printf("Skipping 3 channel calibration.\n");
  333.     else printf("3 channels will be calibrated.\n");
  334.     
  335.     frames=0;
  336.     if(automatic){
  337.         frames=NFRAME;
  338.         printf("Taking photometer readings automatically.\n"
  339.             "Assuming %.3f milliVolts per %s\n",1000.0/LuminancePerVoltage,LP->units);
  340.     }
  341.     else printf("Please enter the photometer readings,in %s, as prompted.\n",LP->units); 
  342.     printf("\n");
  343.  
  344.     printf("First let's check the photometer's gain setting.\n");
  345.     printf("Point the photometer at the screen and hit return:");
  346.     gets(string);
  347.     do{
  348.         SetPort(window);
  349.         BringToFront(window);
  350.         FillRect(&window->portRect,qd.black);
  351.         darkLuminance=GetALuminance(LP,device,frames,LuminancePerVoltage,0,LP->VMax,LP->VMax,LP->VMax);
  352.         FillRect(&window->portRect,qd.white);
  353.         luminance=GetALuminance(LP,device,frames,LuminancePerVoltage,0,LP->VMax,LP->VMax,LP->VMax);
  354.         luminance-=darkLuminance;
  355.         SendBehind(window,NULL);
  356.         printf("White seems to be %.1f %s brighter than black.\n",luminance,LP->units);
  357.         myChar='N';
  358.         printf("Do you want to try that again? (%c):",myChar);
  359.         gets(string);
  360.         sscanf((char*)string,"%c",&myChar);
  361.     }while(toupper(myChar)=='Y');
  362.  
  363.     printf("Please occlude the photometer so that I can read the dark level.\n"
  364.     "Hit cr when ready:");
  365.     gets(string);
  366.     darkLuminance=GetALuminance(LP,device,frames,LuminancePerVoltage,0,LP->VMax,LP->VMax,LP->VMax);
  367.     printf("The dark level of %.2f cd/m^2 will be subtracted from all readings\n"
  368.         ,darkLuminance);
  369.     
  370.     testSize=160;
  371.     printf("Enter, in pixels, the diameter of the test patch (%d):",testSize);
  372.     gets(string);
  373.     sscanf((char *)string,"%d",&testSize);
  374.  
  375.     /* draw circle for placement of photocell */
  376.     SetPort(window);
  377.     myRect = window->portRect;
  378.     InsetRect(&myRect,(myRect.right-testSize)/2,(myRect.bottom-testSize)/2);
  379.     PmForeColor(255);
  380.     PenSize(10,10);
  381.     FrameOval(&myRect);
  382.     SetPort(oldWindow);
  383.     
  384.     printf("Please remove the occluder. \n");
  385.     printf("Place your photometer so as to read the luminance of the "
  386.         "center of the test screen.\n"); 
  387.  
  388.     if(skipRGB)cycles=1;
  389.     else cycles=4;
  390.     printf("Then enter number of times you wish to repeat the whole cycle of measurement. (%d):"
  391.         ,cycles);
  392.     gets((char *)string);
  393.     sscanf((char *)string,"%d",&cycles);
  394.     if(cycles<=0)return;
  395.  
  396.     SetPort(window);
  397.     BringToFront(window);
  398.     FillRect(&window->portRect,qd.gray);
  399.     luminance=GetALuminance(LP,device,frames,LuminancePerVoltage,0,LP->VMax,LP->VMax,LP->VMax);
  400.     luminance-=darkLuminance;
  401.     printf("50%% dithered gray is %.1f %s\n",luminance,LP->units);
  402.     FillRect(&window->portRect,qd.white);
  403.     luminance=GetALuminance(LP,device,frames,LuminancePerVoltage,0,LP->VMax,LP->VMax,LP->VMax);
  404.     luminance-=darkLuminance;
  405.     printf("White is %.1f %s\n",luminance,LP->units);
  406.     FillRect(&window->portRect,qd.black);
  407.     luminance=GetALuminance(LP,device,frames,LuminancePerVoltage,0,LP->VMax,LP->VMax,LP->VMax);
  408.     luminance-=darkLuminance;
  409.     printf("Black is %.1f %s\n",luminance,LP->units);
  410.     SendBehind(window,NULL);
  411.  
  412.     printf("It is essential that the calibration be made with the same background\n"
  413.     "luminance as will ultimately be used in your experiments. You may specify the\n"
  414.     "background luminance either as a luminance in %s, or as a corresponding \n"
  415.     "DAC value (0..%d).\n",LP->units,LP->VMax);
  416.     myChar='Y';
  417.     printf("Do you wish to enter a luminance? (%c):",myChar);
  418.     gets(string);
  419.     sscanf((char *)string,"%c",&myChar);
  420.     if(toupper(myChar)=='Y'){
  421.         printf("Enter the desired background luminance in %s (%f):",LP->units,LP->LBackground);
  422.         gets(string);
  423.         sscanf(string,"%lf",&LP->LBackground);
  424.         printf("Now hunting for the DAC value . . .\n");
  425.         LToVHunt(LP,device,(CWindowPtr)window,LuminancePerVoltage,frames,darkLuminance);
  426.         printf("Will use the nearest DAC value, which is %d\n",LP->VBackground);
  427.     }
  428.     else {
  429.         printf("Enter the DAC value for the background (%d):",LP->VBackground);
  430.         gets(string);
  431.         sscanf((char *)string,"%d",&LP->VBackground);
  432.     }
  433.     LP->table[2].rgb.red = LP->table[2].rgb.green =  LP->table[2].rgb.blue = 
  434.         LP->VBackground<<LP->leftShift;
  435.     LoadLuminances(device,LP,2,2);
  436.     SetPort(window);
  437.     BringToFront(window);
  438.     PmBackColor(2);
  439.     EraseRect(&window->portRect);
  440.     luminance=GetALuminance(LP,device,frames,LuminancePerVoltage,0,LP->VMax,LP->VMax,LP->VMax);
  441.     luminance-=darkLuminance;
  442.     SendBehind(window,NULL);
  443.     SetPort(oldWindow);
  444.     printf("Your background luminance is %.3f %s\n",luminance,LP->units);
  445.  
  446.     /* plan the luminance and gain measurements */
  447.     if(skipRGB)nLrgb=0;
  448.     readTotal=cycles*(ROUND_LUMINANCES+3*nLrgb);
  449.     printf("Now taking %d readings . . .\n",readTotal);
  450.     for (i=0;i<nL;i++) L[i]=0.0;
  451.     k=0;
  452.     for (i=0;i<nL;i++){
  453.         V[i]=k;
  454.         if(k==LP->VMax) {
  455.             if(i+1<nL) V[++i]=k;    /* put two points at VMax, because it's important */
  456.             nL=i+1;
  457.             break;
  458.         }
  459.         k += (LP->VMax+1)/(nL-2);
  460.         if(k>LP->VMax)k=LP->VMax;
  461.     }
  462.     for (k=0;k<3;k++)for(i=0;i<nLrgb;i+=2){
  463.         for(j=0;j<3;j++) nrgb[k][i+1][j]=nrgb[k][i][j]=(LP->VMax+1)*(0.667+0.333*i/nLrgb);
  464.         nrgb[k][i+1][k]=LP->VMax;
  465.         nrgb[k][i][k]=0;
  466.         if(k==2) nrgb[k][i][k]=(LP->VMax+1)/2;    /* take smaller step, because blue gain is highest */
  467.     }
  468.     for(k=0;k<3;k++)for (i=0;i<nLrgb;i++)Lrgb[k][i]=0.0;
  469.     
  470.     /* measure several times to minimize the effect of luminance drift */
  471.     SetPort(window);
  472.     BringToFront(window);
  473.     GetFNum((StringPtr) "\pMonaco",&FontNum);
  474.     TextFont(FontNum);
  475.     TextSize(36);
  476.     myRect = window->portRect;
  477.     bottom=myRect.bottom;
  478.     SetRect(&TestRect, (myRect.right+myRect.left-testSize)/2
  479.         ,(myRect.top+myRect.bottom-testSize)/2
  480.         ,(myRect.right+myRect.left+testSize)/2
  481.         ,(myRect.top+myRect.bottom+testSize)/2 );
  482.     
  483.     /* set up background */
  484.     LP->table[2].rgb.red = LP->table[2].rgb.green =  LP->table[2].rgb.blue = 
  485.         LP->VBackground<<LP->leftShift;
  486.     LoadLuminances(device,LP,1,2);    /* load background into clut */
  487.     PmBackColor(2);                /* our background */
  488.     EraseRect(&myRect);
  489.     PmBackColor(1);                /* our test luminance */
  490.     EraseOval(&TestRect);
  491.     FlushEvents(everyEvent,0);
  492.     readNumber = 0;    /* initialize to count the number of reading */
  493.     for(j=0;j<cycles;j++)    {
  494.         for (k=-1;k<3;k++)    {
  495.             if(k<0) ii=nL;
  496.             else ii=nLrgb;
  497.             for (i=0;i<ii;i++){
  498.                 /* set up test */
  499.                 SetRect(&labelRect,0,bottom-40,300,bottom);
  500.                 SetPort(window);
  501.                 PmBackColor(2);
  502.                 EraseRect(&labelRect);
  503.                 PmForeColor(ctSize);                /* black */
  504.                 MoveTo(0,bottom);
  505.                 if(k<0) DrawPrintf("%3d",V[i]);
  506.                 else DrawPrintf("%3d%4d%4d",k,nrgb[k][i][0],nrgb[k][i][1],nrgb[k][i][2]);
  507.                 SetPort(oldWindow);
  508.                 readNumber++;
  509.                 if(!automatic){
  510.                     printf("%d out of %d readings:\t",readNumber,readTotal);    
  511.                 }
  512.                 if(k<0) {
  513.                     luminance=GetALuminance(LP,device,frames,LuminancePerVoltage,1
  514.                         ,V[i],V[i],V[i]);
  515.                     luminance-=darkLuminance;
  516.                     L[i] += luminance/cycles;
  517.                 }
  518.                 else {
  519.                     luminance=GetALuminance(LP,device,frames,LuminancePerVoltage,1
  520.                         ,nrgb[k][i][0],nrgb[k][i][1],nrgb[k][i][2]);
  521.                     luminance-=darkLuminance;
  522.                     Lrgb[k][i] += luminance/cycles;
  523.                 }
  524.             }
  525.         }
  526.     }
  527.     GDDisposeWindow1(window);
  528.     SetPort(oldWindow);
  529.     SetGDevice(oldDevice);
  530.     
  531.     /* polynomial and quadratic fits */
  532.     ma=MAX_COEFFICIENTS;
  533.     x=vector(1,nL);
  534.     y=vector(1,nL);
  535.     sig=vector(1,nL);
  536.     a=vector(1,nL);
  537.     u=matrix(1,nL,1,ma);
  538.     v=matrix(1,ma,1,ma);
  539.     w=vector(1,ma);
  540.     cvm=matrix(1,ma,1,ma);
  541.     for(i=0;i<nL;i++){
  542.         x[i+1]=V[i];
  543.         y[i+1]=L[i];
  544.         sig[i+1]=10.0;
  545.     }
  546.     svdfit(x,y,sig,nL,a,ma,u,v,w,&chisq,fpoly);    /* ma-1th order polynomial curve fit */
  547.     svdvar(v,ma,w,cvm);
  548.     if(ma>MAX_COEFFICIENTS)PrintfExit("Error: too many coefficients\007\n");
  549.     LP->coefficients=ma;
  550.     for(i=0;i<ma;i++) LP->p[i]=a[i+1];
  551.     printf("L(V) =");
  552.     for(i=0;i<ma;i++) printf(" + %6g V^%d",a[i+1],i);
  553.     printf(".  chisq %g\n",chisq);
  554.     svdfit(x,y,sig,nL,a,3,u,v,w,&chisq,fpoly);    /* 2nd order polynomial curve fit */
  555.     svdvar(v,3,w,cvm);
  556.     printf("\nquadratic fit:\n\n");
  557.     for(i=0;i<3;i++) LP->q[i]=a[i+1];
  558.     printf("L(V) =");
  559.     for(i=0;i<3;i++) printf(" + %g V^%d",a[i+1],i);
  560.     printf(".  chisq %g\n",chisq);
  561.     free_vector(x,1,nL);
  562.     free_vector(y,1,nL);
  563.     free_vector(sig,1,nL);
  564.     free_vector(a,1,nL);
  565.     free_matrix(u,1,nL,1,ma);
  566.     free_matrix(v,1,ma,1,ma);
  567.     free_vector(w,1,ma);
  568.     free_matrix(cvm,1,ma,1,ma);
  569.     e=0.0;
  570.     for(i=0;i<nL;i++){
  571.         f=0.0;
  572.         VV=1.0;
  573.         for(j=0;j<LP->coefficients;j++){
  574.             f+=LP->p[j]*VV;
  575.             VV*=V[i];
  576.         }
  577.         e+=(L[i]-f)*(L[i]-f);
  578.     }
  579.     LP->polynomialError = sqrt(e/nL);
  580.     e=0.0;
  581.     for(i=0;i<nL;i++){
  582.         f=0.0;
  583.         VV=1.0;
  584.         for(j=0;j<3;j++){
  585.             f+=LP->q[j]*VV;
  586.             VV*=V[i];
  587.         }
  588.         e+=(L[i]-f)*(L[i]-f);
  589.     }
  590.     LP->quadraticError = sqrt(e/nL);
  591.     
  592.     /* power law fit */
  593.     /* L=p[1]+Rectify(p[2]+p[3]*V)^p[4] */
  594.     p=vector(1,4);    /* initial starting point */
  595.     /*
  596.         It is necessary to have a reasonable starting point or the search can get
  597.         stuck out in the boondocks. Since the search is quite slow it is desirable to
  598.         give it as good a starting point as possible. Therefore I make use of the
  599.         quadratic fit to try to give it a pretty good starting point.
  600.     */
  601.     /* use the quadratic fit as the starting point */
  602.     p[1]=LP->q[0]-0.25*LP->q[1]*LP->q[1]/LP->q[2];
  603.     p[2]=0.5*LP->q[1]/sqrt(LP->q[2]);
  604.     p[3]=sqrt(LP->q[2]);
  605.     p[4]=2.0;
  606.     /*
  607.         On second thought, try starting with a gamma of 2.3.
  608.         Set the luminance offset p[1] to the measured dc level.
  609.         Set the brightness p[2] so as to match the number offset in the quadratic fit.
  610.         Set the contrast p[3] so as to fit the maximum data point.
  611.         Set gamma p[4] to 2.3
  612.     */
  613.     p[1]=L[1];
  614.     p[2]=0.5*LP->q[1]/sqrt(LP->q[2]);
  615.     p[4]=2.3;
  616.     p[3]=(pow(L[nL-1]-p[1],1.0/p[4])-p[2])/V[nL-1];
  617.     p[2]=p[3]*0.5*LP->q[1]/LP->q[2];
  618.     p[3]=(pow(L[nL-1]-p[1],1.0/p[4])-p[2])/V[nL-1];
  619.     xi=matrix(1,4,1,4);    /* initial set of directions */
  620.     for(i=1;i<=4;i++)for(j=1;j<=4;j++)xi[i][j]=0.0;
  621.     xi[1][1]=0.1;
  622.     xi[2][2]=1.0;
  623.     xi[3][3]=0.1;
  624.     xi[4][4]=0.01;
  625.     ftol=1e-2;    /* fractional tolerance on function value when done */
  626.     printf("\npower law fit:\n\n");
  627.     variables=4;    /* inform PowerRMSError */
  628.     fret=PowerRMSError(p);
  629.     printf("L(V) = %9.4f + Rectify(%9.4f + %9.4f V)^%9.4f ± %9.4f\n",p[1],p[2],p[3],p[4],fret);
  630.     powell(p,xi,4,ftol,&iter,&fret,&PowerRMSError);
  631.     printf("L(V) = %9.4f + Rectify(%9.4f + %9.4f V)^%9.4f ± %9.4f\n",p[1],p[2],p[3],p[4],fret);
  632.     for(i=0;i<4;i++) LP->power[i]=p[i+1];
  633.     free_vector(p,1,4);
  634.     free_matrix(xi,1,4,1,4);
  635.     e=0.0;
  636.     for(i=0;i<nL;i++){
  637.         f=LP->power[1]+LP->power[2]*V[i];
  638.         if(f>0.0) f=LP->power[0]+pow(f,LP->power[3]);
  639.         else f=LP->power[0];
  640.         e+=(L[i]-f)*(L[i]-f);
  641.     }
  642.     LP->powerError = sqrt(e/nL);
  643.     LP->VMin=0;                    /* minimum value that can be loaded into DAC */
  644.     LP->VMax=(1L<<LP->dacSize)-1;    /* maximum value that can be loaded into DAC */
  645.     LP->LMin=VToL(LP,LP->VMin);    /* min luminance */
  646.     LP->LMax=VToL(LP,LP->VMax);    /* max luminance */
  647.     
  648.     /* power law fit with a FIXED gamma */
  649.     /* This is solely for study of the effects of the contrast and brightness knobs */
  650.     /* L=p[1]+Rectify(p[2]+p[3]*V)^p[4] */
  651.     p=vector(1,4);    /* initial starting point */
  652.     /*
  653.         It is necessary to have a reasonable starting point or the search can get
  654.         stuck out in the boondocks. Since the search is quite slow it is desirable to
  655.         give it as good a starting point as possible. 
  656.     */
  657.     for(i=0;i<4;i++)p[i+1]=LP->power[i];
  658.     p[4]=2.28;
  659.     xi=matrix(1,3,1,3);    /* initial set of directions */
  660.     for(i=1;i<=3;i++)for(j=1;j<=3;j++)xi[i][j]=0.0;
  661.     xi[1][1]=0.1;
  662.     xi[2][2]=1.0;
  663.     xi[3][3]=0.1;
  664.     ftol=1e-2;        /* fractional tolerance on function value when done */
  665.     variables=3;    /* global, to inform PowerRMSError() */
  666.     printf("\npower law fit, with fixed gamma:\n\n");
  667.     fret=PowerRMSError(p);
  668.     printf("L(V) = %9.4f + Rectify(%9.4f + %9.4f V)^%9.4f ± %9.4f\n",p[1],p[2],p[3],p[4],fret);
  669.     powell(p,xi,3,ftol,&iter,&fret,&PowerRMSError);
  670.     printf("L(V) = %9.4f + Rectify(%9.4f + %9.4f V)^%9.4f ± %9.4f\n",p[1],p[2],p[3],p[4],fret);
  671.     for(i=0;i<4;i++) LP->fixedPower[i]=p[i+1];
  672.     free_vector(p,1,4);
  673.     free_matrix(xi,1,3,1,3);
  674.     e=0.0;
  675.     for(i=0;i<nL;i++){
  676.         f=LP->fixedPower[1]+LP->fixedPower[2]*V[i];
  677.         if(f>0.0) f=LP->fixedPower[0]+pow(f,LP->fixedPower[3]);
  678.         else f=LP->fixedPower[0];
  679.         e += (L[i]-f)*(L[i]-f);
  680.     }
  681.     LP->fixedPowerError = sqrt(e/nL);
  682.     
  683.     if(skipRGB){
  684.         printf("Please enter gain of red channel   (%6.4f):",LP->r);
  685.         gets(string);
  686.         sscanf((char *)string,"%lf",&LP->r);
  687.         printf("Please enter gain of green channel (%6.4f):",LP->g);
  688.         gets(string);
  689.         sscanf((char *)string,"%lf",&LP->g);
  690.         printf("Please enter gain of blue channel  (%6.4f):",LP->b);
  691.         gets(string);
  692.         sscanf((char *)string,"%lf",&LP->b);
  693.         gain[0]=LP->r;
  694.         gain[1]=LP->g;
  695.         gain[2]=LP->b;
  696.     }
  697.     else {    
  698.         /* calculate R,G,B gains */
  699.         for(k=0;k<3;k++){
  700.             gain[k]=0.0;
  701.             for(i=0;i<nLrgb;i+=2)
  702.                 gain[k]+=(LToV(LP,Lrgb[k][i])-LToV(LP,Lrgb[k][i+1]))/(nrgb[k][i][k]-nrgb[k][i+1][k]);
  703.             gain[k] /= nLrgb/2;
  704.             printf("gain[%d]=%6.4f\n",k,gain[k]);
  705.         }
  706.     }
  707.     printf("Sum of gains %6.4f\n",gain[0]+gain[1]+gain[2]);
  708.     LP->r=gain[0]/(gain[0]+gain[1]+gain[2]);        /* Normalize gain */
  709.     LP->g=gain[1]/(gain[0]+gain[1]+gain[2]);
  710.     LP->b=gain[2]/(gain[0]+gain[1]+gain[2]);
  711.     e=1.0-(gain[0]+gain[1]+gain[2]);
  712.     if(skipRGB && fabs(e)<1e-6){
  713.         printf("Please enter gain accuracy   (%6.4f):",LP->gainAccuracy);
  714.         gets(string);
  715.         sscanf((char *)string,"%lf",&LP->gainAccuracy);
  716.     }
  717.     else LP->gainAccuracy=e;
  718.     LP->LBackground=VToL(LP,LP->VBackground);
  719.     LP->gm=(LP->VMax/(2.0*LP->LBackground))
  720.         *(LP->LBackground-VToL(LP,LToV(LP,LP->LBackground)-1.0));
  721.     
  722.     SaveData(LP,nL,V,L,nLrgb,nrgb,Lrgb,gain);
  723.     SaveLuminanceRecord(LP);
  724.     printf(
  725.         "Congratulations! You've finished the luminance calibration of screen %d. The calibration results\n"
  726.         "have been saved as CalibrateLuminance%d.data and LuminanceRecord%d.h.\n\n"
  727.         
  728.         "The CalibrateLuminance%d.data file is a text file suitable for graphing the monitor's gamma function.\n"
  729.         "We recommend using CricketGraph with the CalibrateLuminance.format. The V column is the equivalent\n"
  730.         "number loaded into the clut (V=%d*v), and the L column is the measured luminance.\n\n"
  731.         
  732.         "The program's results have also been stored in a C header file called LuminanceRecord%d.h that all\n"
  733.         "your programs should read by using either #include at compile time, or ReadLuminanceRecord() at\n"
  734.         "runtime. LuminanceRecord%d.h describes the gains of the three video pathways of your video card and\n"
  735.         "ISR Video Attenuator, and describes the gamma function of your monitor, to allow automatic gamma\n"
  736.         "correction later, using SetLuminance(), etc. See Pelli & Zhang (1991) Vision Research 31:1337-1360.\n"
  737.         "Hit RETURN to exit. Good luck.\n"
  738.         ,LP->screen,LP->screen,LP->screen,LP->screen,LP->VMax,LP->screen,LP->screen
  739.     );
  740. }
  741.  
  742. void SaveData(luminanceRecord *LP,int nL,int V[],double L[]
  743.     ,int nLrgb,int nrgb[3][NLRGB][3],double Lrgb[3][NLRGB],double gain[3])
  744. /* create Cricket data file, to graph the gamma function */
  745. {
  746.     FILE *myfile;
  747.     int i,j,k;
  748.     double VV,LFit;
  749.     static char fileName[80];
  750.     
  751.     sprintf(fileName,"LuminanceCalibration%d.data",LP->screen);
  752.     myfile=fopen(fileName,"w");
  753.     SetFileInfo(fileName,'TEXT','CGRF');
  754.     fprintf(myfile,"*\n");
  755.     fprintf(myfile,"V\tL \tV(L) ");
  756.     fprintf(myfile,"\tquadratic\tpolynomial\tpower ");
  757.     fprintf(myfile,"\tL-quadratic\tL-polynomial\tL-power ");
  758.     fprintf(myfile,"\tVr \tVg \tVb ");
  759.     fprintf(myfile,"\tLr \tLg \tLb \tV(Lr) \tV(Lg) \tV(Lb) ");
  760.     fprintf(myfile,"\tLR.r\tLR.g \tLR.b \tsum of gains");
  761.     for(i=0;i<nL;i++)    {
  762.         fprintf(myfile,"\n%d\t",V[i]);
  763.         fprintf(myfile,"%9.4g\t",L[i]);
  764.         fprintf(myfile,"%9.4g\t",LToV(LP,L[i]));
  765.         /* quadratic fit */
  766.         LFit=0.0;
  767.         VV=1.0;
  768.         for(j=0;j<3;j++){
  769.             LFit += LP->q[j]*VV;
  770.             VV *= V[i];
  771.         }
  772.         fprintf(myfile,"%9.4g\t",LFit);
  773.         /* polynomial fit */
  774.         LFit=0.0;
  775.         VV=1.0;
  776.         for(j=0;j<LP->coefficients;j++){
  777.             LFit += LP->p[j]*VV;
  778.             VV *= V[i];
  779.         }
  780.         fprintf(myfile,"%9.4g\t",LFit);
  781.         /* power law fit */
  782.         LFit=LP->power[1]+LP->power[2]*V[i];
  783.         if(LFit>0.0) LFit=LP->power[0]+pow(LFit,LP->power[3]);
  784.         else LFit=LP->power[0];
  785.         fprintf(myfile,"%9.4g\t",LFit);
  786.  
  787.         /* quadratic fit */
  788.         LFit=0.0;
  789.         VV=1.0;
  790.         for(j=0;j<3;j++){
  791.             LFit += LP->q[j]*VV;
  792.             VV *= V[i];
  793.         }
  794.         fprintf(myfile,"%9.4g\t",L[i]-LFit);
  795.         /* polynomial fit */
  796.         LFit=0.0;
  797.         VV=1.0;
  798.         for(j=0;j<LP->coefficients;j++){
  799.             LFit += LP->p[j]*VV;
  800.             VV *= V[i];
  801.         }
  802.         fprintf(myfile,"%9.4g\t",L[i]-LFit);
  803.         /* power law fit */
  804.         LFit=LP->power[1]+LP->power[2]*V[i];
  805.         if(LFit>0.0) LFit=LP->power[0]+pow(LFit,LP->power[3]);
  806.         else LFit=LP->power[0];
  807.         fprintf(myfile,"%9.4g\t",L[i]-LFit);
  808.  
  809.         if(i<nLrgb) for(k=0;k<3;k++) {
  810.             VV=LP->r*nrgb[k][i][0]+LP->g*nrgb[k][i][1]+LP->b*nrgb[k][i][2];
  811.             fprintf(myfile,"%f\t",VV);
  812.         }
  813.         if(i<nLrgb) for(k=0;k<3;k++) fprintf(myfile,"%f\t",Lrgb[k][i]);
  814.         if(i<nLrgb) for(k=0;k<3;k++) fprintf(myfile,"%f\t",LToV(LP,Lrgb[k][i]));
  815.         if(i==0) fprintf(myfile,"%6g\t%6g\t%6g\t%6g",LP->r,LP->g,LP->b
  816.                 ,gain[0]+gain[1]+gain[2]);
  817.     }
  818.     fclose(myfile);
  819. }
  820.  
  821. void SaveLuminanceRecord(luminanceRecord *LP)
  822. /* Create C include file "LuminanceRecord?.h", where ? is the screen number. */
  823. {
  824.     FILE *file;
  825.     int i;
  826.     static char string[100],dateString[100];
  827.     long seconds;
  828.     
  829.     LP->VMin=0;
  830.     LP->VMax=(1L<<LP->dacSize)-1;    /* maximum value that can be loaded into DAC */
  831.     LP->LMin=VToL(LP,LP->VMin);        /* min luminance */
  832.     LP->LMax=VToL(LP,LP->VMax);        /* max luminance */
  833.     sprintf(string,"LuminanceRecord%d.h",LP->screen);
  834.     GetDateTime((void *)&seconds);
  835.     file=fopen(string,"w");
  836.     SetFileInfo(string,'TEXT','KAHL');
  837.     printf("Please enter the following descriptive information for the %s file.\n",string);
  838.  
  839.     fprintf(file,"/* This %s file is a description of a monitor produced by the CalibrateLuminance program. */\n",string);
  840.     fprintf(file,"/* This file may be #included in any C program. It simply fills in a luminanceRecord data structure. */\n");
  841.     fprintf(file,"/* Or use ReadLuminanceRecord.c to read this file at runtime. */\n");
  842.     fprintf(file,"/* The theory is described by Pelli & Zhang (1991). The data structure is defined in Luminance.h.  */\n");
  843.     fprintf(file,"/* CalibrateLuminance and Luminance.h are part of the VideoToolbox software. */\n");
  844.     fprintf(file,"/* Pelli & Zhang (1991) Accurate control of contrast on microcomputer displays. */\n");
  845.     fprintf(file,"/* Vision Research, 31:1337-1360. */\n");
  846.     fprintf(file,"/* The VideoToolbox software is available, free for research purposes, from Denis Pelli. */\n");
  847.     fprintf(file,"/* Institute for Sensory Research, Syracuse University, Syracuse, NY 13244-5290. */\n");
  848.     fprintf(file,"LR.screen=%d;    /* device=GetScreenDevice(LR.screen); */\n"
  849.         ,LP->screen);
  850.     fprintf(file,"/* Caution: the screen number used here and in GetScreen Device is NOT the same as */\n");
  851.     fprintf(file,"/* displayed by the Monitors cdev in the Control Panel. Sorry. The most obvious difference */\n");
  852.     fprintf(file,"/* is that GetScreenDevice always assigns 0 to the main screen, the one with the menu bar. */\n");
  853.     /* time and date */
  854.     IUDateString(seconds,longDate,(unsigned char *)string);
  855.     PtoCstr((unsigned char *)string);
  856.     IUTimeString(seconds,FALSE,(unsigned char *)dateString);
  857.     PtoCstr((unsigned char *)dateString);
  858.     sprintf(dateString,"%s %s",dateString,string);
  859.     LP->date=dateString;
  860.     fprintf(file,"LR.date=\"%s\";\n",LP->date);
  861.     printf("LR.date=\"%s\";\n",LP->date);    
  862.  
  863.     printf("Enter LR.id, i.e. monitor model and serial # (%s):",LP->id);
  864.     gets(string);
  865.     if(strlen(string)==0) strcpy(string,LP->id);
  866.     printf("%s\n",string);
  867.     fprintf(file,"LR.id=\"%s\";\n",string);
  868.  
  869.     printf("Enter LR.name, i.e. informal monitor name \n(%s):",LP->name);
  870.     gets(string);
  871.     if(strlen(string)==0) strcpy(string,LP->name);
  872.     printf("%s\n",string);
  873.     fprintf(file,"LR.name=\"%s\";\n",string);
  874.  
  875.     printf("Enter LR.notes, i.e. details of calibration: who & how\n(%s):\n",LP->notes);
  876.     gets(string);
  877.     if(strlen(string)==0) strcpy(string,LP->notes);
  878.     printf("%s\n",string);
  879.     fprintf(file,"LR.notes=\"%s\";\n",string);
  880.  
  881.     printf("Nominally, the Apple High-Res. Monochrome Monitor has 76 pixels/inch.\n");
  882.     printf("Nominally, the Apple High-Res. RGB Monitor has 69 pixels/inch.\n");
  883.     printf("Enter LR.dpi, i.e. pixels per inch (%.1f):",LP->dpi);
  884.     gets(string);
  885.     sscanf(string,"%lf",&LP->dpi);
  886.     printf("%.1f\n",LP->dpi);
  887.     fprintf(file,"LR.dpi=%.1f;\t/* pixels per inch */\n",LP->dpi);
  888.  
  889.     fprintf(file,"LR.Hz=%.2f;\t/* frames per second */\n",LP->Hz);
  890.  
  891.     printf("Enter LR.units luminance units (%s):",LP->units);
  892.     gets(string);
  893.     if(strlen(string)==0) strcpy(string,LP->units);
  894.     printf("%s\n",string);
  895.     fprintf(file,"LR.units=\"%s\";\n",string);
  896.  
  897.     fprintf(file,"/* coefficients of polynomial fit */\n");
  898.     fprintf(file,"LR.coefficients=%ld;    /* # of coefficients in polynomial fit */\n"
  899.         ,LP->coefficients);
  900.     fprintf(file,"/* L(V)=p[0]+p[1]*V+p[2]*V*V+ . . . ±polynomialError */\n");
  901.     for(i=0;i<LP->coefficients;i++)    {
  902.         fprintf(file,"LR.p[%d]=%6g;\n",i,LP->p[i]);
  903.     }
  904.     fprintf(file,"LR.polynomialError=%8.4f;    /* RMS error of fit */\n",LP->polynomialError);
  905.     fprintf(file,"/* coefficients of quadratic fit */\n");
  906.     fprintf(file,"/* L(V)=q[0]+q[1]*V+q[2]*V*V±quadraticError */\n");
  907.     for(i=0;i<3;i++)    {
  908.         fprintf(file,"LR.q[%d]=%6g;\n",i,LP->q[i]);
  909.     }
  910.     fprintf(file,"LR.quadraticError=%8.4f;    /* RMS error of fit */\n",LP->quadraticError);
  911.     fprintf(file,"/* coefficients of power law fit */\n");
  912.     fprintf(file,"/* L(V)=power[0]+Rectify(power[1]+power[2]*V)^power[3]±powerError */\n");
  913.     fprintf(file,"/* where Rectify(x)=x if x≥0, and Rectify(x)=0 otherwise */\n");
  914.     fprintf(file,"/* Pelli & Zhang (1991) Eqs.9&10 use symbols v=V/VMax, alpha=power[0], beta=power[1], kappa=power[2]*255, gamma=power[3] */\n");
  915.     for(i=0;i<4;i++)    {
  916.         fprintf(file,"LR.power[%d]=%6g;\n",i,LP->power[i]);
  917.     }
  918.     fprintf(file,"LR.powerError=%8.4f;    /* RMS error of fit */\n",LP->powerError);
  919.     fprintf(file,"/* coefficients of power law fit, with fixed exponent */\n");
  920.     fprintf(file,"/* L(V)=fixedPower[0]+Rectify(fixedPower[1]+fixedPower[2]*V)^fixedPower[3]±fixedPowerError */\n");
  921.     for(i=0;i<4;i++)    {
  922.         fprintf(file,"LR.fixedPower[%d]=%6g;\n",i,LP->fixedPower[i]);
  923.     }
  924.     fprintf(file,"LR.fixedPowerError=%8.4f;    /* RMS error of fit */\n",LP->fixedPowerError);
  925.     fprintf(file,"LR.r=%6g;\n",LP->r);
  926.     fprintf(file,"LR.g=%6g;\n",LP->g);
  927.     fprintf(file,"LR.b=%6g;\n",LP->b);
  928.     fprintf(file,"LR.gainAccuracy=%6g;\n",LP->gainAccuracy);
  929.     fprintf(file,"LR.gm=%6g;    /* The monitor's contrast gain. */\n",LP->gm);
  930.     fprintf(file,"LR.dacSize=%d;    /* bits */\n",LP->dacSize);
  931.     fprintf(file,"LR.VMin=%3d;    /* minimum value that can be loaded into DAC */\n",LP->VMin);
  932.     fprintf(file,"LR.VMax=%3d;    /* maximum value that can be loaded into DAC */\n",LP->VMax);
  933.     fprintf(file,"LR.LMin=%8.2f;    /* luminance at VMin */\n",LP->LMin);
  934.     fprintf(file,"LR.LMax=%8.2f;    /* luminance at VMax */\n",LP->LMax);
  935.     fprintf(file,"LR.LBackground=%8.3f;    /* background luminance during calibration */\n",LP->LBackground);
  936.     fprintf(file,"LR.VBackground=%d;    /* background number used during calibration */\n",LP->VBackground);
  937.     fprintf(file,"LR.rangeSet=0;    /* indicate that range parameters have yet to be set */\n");
  938.     fprintf(file,"LR.L.exists=0;    /* indicate that luminance table has yet to be initialized */\n");
  939.     fclose(file);
  940. }
  941.  
  942. FLOAT PowerRMSError(FLOAT pp[])
  943. /*
  944. Returns average squared error of power law fit.
  945. pp[] contains the parameters of a power law fit to the
  946. luminance measurements L[] made at numbers V[].
  947. */
  948. {
  949.     int i;
  950.     FLOAT f,e,a;
  951.     extern double L[];
  952.     extern int V[];
  953.     extern int nL;
  954.     extern int variables;
  955.     extern FLOAT *p;
  956.     static FLOAT q[5];
  957.     
  958.     for(i=1;i<=variables;i++)q[i]=pp[i];    /* copy the variable parameters */
  959.     for(i=variables+1;i<=4;i++)q[i]=p[i];    /* copy the fixed parameters */
  960.     e=0.0;
  961.     for(i=0;i<nL;i++){
  962.         f=q[2]+q[3]*V[i];
  963.         if(f>0.0) f=q[1]+pow(f,q[4]);
  964.         else f=q[1];
  965.         a = f-L[i];
  966.         e += a*a;
  967.     }
  968.     e=sqrt(e/nL);
  969.     return e;
  970. }
  971.  
  972.